home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / hist_2d.pro < prev    next >
Text File  |  1997-07-08  |  4KB  |  121 lines
  1. ; $Id: hist_2d.pro,v 1.3 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1992-1997, Research Systems, Inc.  All rights reserved.
  4. ;    Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;    HIST_2D
  8. ;
  9. ; PURPOSE:
  10. ;    Return the density function (histogram) of two variables.
  11. ;
  12. ; CATEGORY:
  13. ;    Image processing, statistics, probability.
  14. ;
  15. ; CALLING SEQUENCE:
  16. ;    Result = hist_2d(V1, V2)
  17. ; INPUTS:
  18. ;    V1 and V2 = arrays containing the variables.  May be any non-complex
  19. ;        numeric type.
  20. ;
  21. ; Keyword Inputs:
  22. ;       MIN1:   MIN1 is the minimum V1 value to consider. If this
  23. ;               keyword is not specified, then it is set to 0.
  24. ;
  25. ;       MIN2:   MIN2 is the minimum V2 value to consider. If this
  26. ;               keyword is not specified, then it is set to 0.
  27. ;
  28. ;       MAX1:   MAX1 is the maximum V1 value to consider. If this
  29. ;               keyword is not specified, then V1 is searched for
  30. ;               its largest value.
  31. ;
  32. ;       MAX2    MAX2 is the maximum V2 value to consider. If this
  33. ;               keyword is not specified, then V2 is searched for
  34. ;               its largest value.
  35. ;
  36. ;       BIN1    The size of each bin in the V1 direction (column
  37. ;               width).  If this keyword is not specified, the
  38. ;               size is set to 1.
  39. ;
  40. ;       BIN2    The size of each bin in the V2 direction (row
  41. ;               height).  If this keyword is not specified, the
  42. ;               size is set to 1.
  43. ;
  44. ; OUTPUTS:
  45. ;    The two dimensional density function of the two variables,
  46. ;    a longword array of dimensions (m1, m2), where:
  47. ;        m1 = Floor((max1-min1)/bin1) + 1
  48. ;       and  m2 = Floor((max2-min2)/bin2) + 1
  49. ;    and Result(i,j) is equal to the number of sumultaneous occurences
  50. ;    of an element of V1 falling in the ith bin, with the same element
  51. ;    of V2 falling in the jth bin, where:
  52. ;        i = (v1 < max1 - min1 > 0) / b1
  53. ;       and  j = (v2 < max2 - min2 > 0) / b2
  54. ;
  55. ;    Note: elements larger than the max or smaller than the min are
  56. ;    truncated to the max and min, respectively.
  57. ;    
  58. ; COMMON BLOCKS:
  59. ;    None.
  60. ; SIDE EFFECTS:
  61. ;    None.
  62. ; RESTRICTIONS:
  63. ;    Not usable with complex or string data.
  64. ; PROCEDURE:
  65. ;    Creates a combines array from the two variables, equal to the
  66. ;    linear subscript in the resulting 2D histogram, then applies
  67. ;    the standard histogram function.
  68. ;
  69. ; EXAMPLE:
  70. ;    Return the 2D histogram of two byte images:
  71. ;        R = HIST_2D(image1, image2)
  72. ;
  73. ;    Return the 2D histogram made from two floating point images
  74. ;    with range of -1 to +1, and with 101 (= 2/.02 + 1) bins:
  75. ;        R = HIST_2D(f1, f2, MIN1=-1, MIN2=-1, MAX1=1, MAX2=1, $
  76. ;            BIN1=.02, BIN2=.02)
  77. ;
  78. ; MODIFICATION HISTORY:
  79. ;     Written by:
  80. ;    DMS, Sept, 1992        Written
  81. ;    DMS, Oct, 1995        Added MIN, MAX, BIN keywords following
  82. ;                suggestion of Kevin Trupie, GSC, NASA/GSFC.
  83. ;-
  84. ; Form the 2 dimensional histogram of two arrays.
  85. ; Result(i,j) = density of value i in im1, and value j in im2.
  86. ; Input images must be, of course, the same size....
  87. ;
  88. function hist_2d, im1, im2, Min1 = mn1, Min2 = mn2, Max1 = mx1, Max2 = mx2, $
  89.     Bin1 = b1, Bin2 = b2
  90.  
  91. m1 = max(im1, min=mm1)    ;Find extents of arrays.
  92. m2 = max(im2, min=mm2)
  93.  
  94. if N_elements(mn1) eq 0 then mn1 = 0    ;Supply default values for keywords.
  95. if N_elements(mx1) eq 0 then mx1 = m1
  96. if N_elements(mn2) eq 0 then mn2 = 0
  97. if N_elements(mx2) eq 0 then mx2 = m2
  98. if N_elements(b1) le 0 then b1 = 1L
  99. if N_elements(b2) le 0 then b2 = 1L
  100.  
  101. m1 = floor((mx1-mn1) / b1) + 1L        ;Get # of bins for each
  102. m2 = floor((mx2-mn2) / b2) + 1L
  103.  
  104. if m1 le 0 or m2 le 0 then message,'Illegal bin size'
  105.  
  106.     ; Combine with im1 in low part & im2 in high
  107. if      mn1 eq 0 and mn2 eq 0 and $    ; Fast case without scaling?
  108.     b1 eq 1L and b2 eq 1L and $
  109.     mx1 le m1 and mx2 le m2 and $
  110.     mm1 ge 0 and mm2 ge 0 then $
  111.    h = m1 * long(im2) + long(im1) $
  112. else if b1 eq 1L and b2 eq 1L then $    ;Avoid dividing?
  113.    h = m1 * long((im2 < mx2) - mn2 > 0L) + long((im1 < mx1) - mn1 > 0L) $
  114. else $                    ;Slowest case
  115.    h = m1 * long(((im2 < mx2) - mn2 > 0L) / b2) + $
  116.         long(((im1 < mx1) - mn1 > 0L) / b1)
  117.  
  118. h = histogram(h, min = 0, max= m1 * m2 -1)  ;Get the 1D histogram
  119. return, reform(h, m1, m2, /overwrite) ;and make it 2D
  120. end
  121.